COMPUTING GROUND STATES OF SPIN-1 BOSE-EINSTEIN 
CONDENSATES BY THE NORMALIZED GRADIENT FLOW 
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Abstract. In this paper, we propose an efficient and accurate numerical method for computing 
the ground state of spin-1 Bose-Einstein condensates (BEG) by using the normalized gradient flow or 
imaginary time method. The key idea is to find a third projection or normalization condition based 
on the relation between the chemical potentials so that the three projection parameters used in the 
projection step of the normalized gradient fiow are uniquely determined by this condition as well as 
the other two physical conditions given by the conservation of total mass and total magnetization. 
This allows us to successfully extend the most popular and powerful normalized gradient flow or 
imaginary time method for computing the ground state of single component BEG to compute the 
ground state of spin-1 BEG. An efficient and accurate discretization scheme, the backward- forward 
Euler sine-pseudospectral method (BFSP), is proposed to discretize the normalized gradient flow. 
Extensive numerical results on ground states of spin-1 BEG with ferromagnetic/antiferromagnetic in- 
teraction and harmonic/optical lattice potential in one/three dimensions are reported to demonstrate 
the efficiency of our new numerical method. 
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1. Introduction. Research in low temperature dilute atomic quantum gases 
remains active for more than ten years after the experimental realizations of Bose- 
Einstein condensation (BEC) in alkali atomic gases in 1995 [51[T21[I7]. Extensive the- 
oretical and experimental studies have been carried out to investigate various novel 
phenomena of the condensates. In earlier BEC experiments, the atoms were confined 
in magnetic trap [21 [El [17], in which the spin degrees of freedom is frozen. The 
particles are described by a scalar model and the wavefunction of the particles is gov- 
erned by the Gross-Pitaevskii equation (GPE) within the mean-field approximation 
[TBI [11] [23] . In recent years, experimental achievement of spin-1 and spin-2 conden- 
sates [lOl [181 HU [Ml [28] offers new regimes to study various quantum phenomena 
that are generally absent in a single component condensate. The spinor condensate 
is achieved experimentally when an optical trap, instead of a magnetic trap, is used 
to provide equal confinement for all hyperfine states. 

The theoretical studies of spinor condensate have been carried out in several 
papers since the achievement of it in experiments [TH [20l [221 121] • In contrast to single 
component condensate, a spin-_F {F £ N) condensate is described by a generalized 
coupled GPEs which consists of 2F + 1 equations, each governing one of the 2F + 1 
hyperfine states (mp — —F, ^F+l, F— 1, F) within the mean-field approximation. 
For a spin-1 condensate, at temperature much lower than the critical temperature Tc, 
the three-components wavefunction \E'(x,t) — {ipi{x,t),ipo{x,t),tp^i{x,t)'^ are well 
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described by the following coupled GPEs EH Hi 133 E] , 

2m 

+C2 ij-i^l, 



ihdtiii{y.,t) 
(1.1) 

ihdtipo{-x,t) 
(1.2) 

iHdttp-iix^t) 
(1.3) 



-V^ + y(x) + (co + C2) (1^1 12 + IV^ol') + (co - C2)|V-l|' 



^2^^^' + ^(x) + (co + C2) (l^iP + l^-iP) + col^ol' 
-2c2 V'-i V'o V'l, 

2m 



V2 + l/(x) + (co + C2) (l^A-lP + I^AoP) + (co - C2Ml\' 



1p-l 



+C2 Ipl Ipl 



Here, x = {x, y, z)^ is the Cartesian coordinate vector, t is time, H is the Planck 
constant, m is the atomic mass, and V{'x) is the external trapping potential. When 
a harmonic trap potential is considered. 



(1.4) 



-irf \ ^1 2 2 I 2 2 I 2 2\ 



with uj^, ujy and uj^ being the trap frequencies in the x-, y- and z-direction, respec- 
tively. / and Re(/) denote the conjugate and real part of the function /, respectively. 
There are two atomic collision terms, co = ^^(oo + 2a2) and C2 = ^§^{0.2 — o-o), ex- 
pressed in terms of the s-wave scattering lengths, oq and 02, for scattering channel of 
total hyperfine spin (anti-parallel spin collision) and spin 2 (parallel spin collision), 
respectively. The usual mean- field interaction, co, is positive for repulsive interaction 
and negative for attractive interaction. The spin-exchange interaction, C2, is posi- 
tive for antiferromagnetic interaction and negative for ferromagnetic interaction. The 
wave function is normalized according to 



(1.5) 



*|p:= [ |*(x,i)pdx= / \M^,t)\U^:^ M'-^' 



l = -l 



l=-l 



where N is the total number of particles in the condensate. 

By introducing the dimensionless variables: i — > t/t^m with = iBhi{ujx, Ljy, ujz}, 

X ^ X fls with as = J — ipi — > ^fNipi/ aV^ (l = —1, 0, 1), we get the dimensionless 



coupled GPEs from ([TII|)-(ir3l) as 
1 



(1.6) 

(1.7) 

idtip-i{-K, t) 
(1.8) 



y(x) + (/3„ + Ps) (IV'iP + IV'oP) + (/?« - PsM- 



+(3s'fp-iipo, 
1_ 



V^(x) + iPn+Ps) (IV-ll' + IV'-lP) + /3„|Vo| 



^0 



1 



-V2 + l/(x) + (/3„ + Ps) (l^-iP + iV-oP) + (/3„ - f3s)\M 

+/3s V'o V'l; 



V'- 



where f3n = ^ ^ Ps - ^ - and F(x) = 1(7^.2 



+ Iz^^) with 7^ = ^, 7y 



and 7^ 



Similar as those in single 
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component BEC [241 [51 [511^. in a disk-shaped condensation, i.e. lo^ ~ ujy and cu^ 3> uJx 
7^: = 1, 7y ~ 1 and 72 ^ 1 with aj„, — cj^), the 3D coupled GPEs ([1.6p - (ll.8p can 
be reduced to a 2D coupled GPEs; and in a cigar-shaped condensation, i.e. ujy ^ 
and ujz ^ (^^=^ 7a; = 1, 7a ^ 1 a-nd 7z ^ 1 with w^, = lu^), the 3D coupled 
GPEs ([1.6p - ([1.8p can be reduced to a ID coupled GPEs. Thus here we consider the 
dimensionless coupled GPEs in d-dimensions {d — 1,2,3): 



idtipiix, t) 
(1.9) 

(1.10) 

«9fV'-i(x, t) 
(1.11) 



-iv^ + y(x) + (/3„ + /3s) (IViP + IV'oP) + iPu - PsM-i\' 



+I3s tp^iipo, 
1_ 



y(x) (/3„ + Ps) (iV-iP + iV'-il') + AJ^ol 



-l-2/3s '0-1 V'o V'l, 

-iv^ + V^(X) + (/?„ + A) (iV'-lP + IVoP) + iPn - PsMlf 

-/3s Vo V'l- 



In the equations above, V{x) is a real-valued potential whose shape is determined 
by the type of system under investigation, /?„ cx iV and /3s cx iV correspond to the 
dimensionless mean-field (spin-independent) and spin-exchange interaction, respec- 
tively. Three important invariants of ([1.9p - ([l.lip are the mass (or normalization) of 
the wave function 
(1.12) 

N{^{;t)):=\\'f{;t)f:^f ^ |Vi(x,i)pdx^iV(*(.,0)) = l, t>0, 

l=-l 

the magnetization (with ~1 < M < 1) 

(1.13) M{^i;t)):^ f [|V;i(x,t)|2-|V;„i(x,<)|2]dx = Af(v|/(.,0))=M 

and the energy per particle 

E{^{;t)) = E (^ivv'zp + v{^m'^ + (/?„ - /3s)ivinv'-ip 

|0ir-K|^-i|^ + 2|^o|'(|^i|' + |^-iP 



, Pn , , ,4 , Pn + Ps 

+/3s (V^-iV-oV^i + V'-i'/^oV'i) 1^ rfx = £;(*(•, 0)), i > 



(1.14) 



A fundamental problem in studying BEC is to find the condensate stationary 
states 3>(x), in particular the ground state which is the lowest energy stationary 
state. In other words, the ground state, <I'g(x), is obtained from the minimization of 
the energy functional subject to the conservation of total mass and magnetization: 
Find ($g e S) such that 



(1.15) 



Eg :=£;($,,) -min £;($), 
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where the nonconvex set S is defined as 
(1.16) 

5*= <!$ = 



r I = 1, / [\M^r - i'^-i(x)p] = M, i?($) < ^1 

JR'' ) 



This is a nonconvex minimization problem. When > and /3„ > |/3s| and 
lim|x|^oo ^(x) = oo, the existence of a minimizer of the nonconvex minimization 
problem (ll.isp follows from the standard theory [25j . For understanding the unique- 
ness question note that E{a ■ cf>g) = -E($g) for all a = (e'^S e*^", e'^-i) with 
6'i + 6'_i = 20q. Thus additional constraints have to be introduced to show the 
uniqueness. 

As derived in [S], by defining the Lagrangian 
(1.17) 

£($, Ai, A) E{^) - {U.r + ||0o|P + U-if - l) - A {Uif - - M) , 



we get the Euler-Lagrange equations associated to the minimization problem p.lSp : 

(m + A) 0i(x) = -^V2 + V{x) + (/3„ + ps) (|0iP + l^oP) + {I3n - f3sM-l\' 
(1.18) :=i/i0i, 



(1.19) 



fl M^) = -2^' + + + (I0l|' + l^-lP) + finlM' 

+2Ps 0_i 00 01 := Hq 00, 



(M - A) 0-1 (x) 

(1.20) 



--V^ + y(x) + (/3„ + /3,) (|0_i|2 + |0o|2) + (/3„ - /3,)|0i|2 
-/3s 00 01 := i^-l 0-1. 



Here /i and A are the Lagrange multipliers (or chemical potentials) of the coupled 
GPEs (|1.9p - (|l.lip . In addition, (|1.18|) - (|1.20|) is also a nonhnear eigenvalue problem 
with two constraints 



(1 



(1.22) 



21) / |ci>(x)p dx = / ^ |0,(x)p dx ^ I 

/ [|0i(x)|2-|0_i(x)p]dx = A/. 



-if • = 



In fact, the nonlinear eigenvalue problem (|1.18[) - (|1.20[) can also be obtained from the 
coupled GPEs (fL9)l - (frTT|) by plugging M^^t) = e-'^'*0((x) (/ = 1,0,-1) with 



(1.23) 



fii + /i-i = 2/10. 



Thus it is also called as time-independent coupled GPEs. In physics literatures, any 
eigenfucntion $ of the nonlinear eigenvalue problem (|1.18p - (|1.20[) under constraints 
()1.2ip and (|1.22p . whose energy is larger than the energy of the ground state is called 
as an excited state of the coupled GPEs (|1.9p - (ll.lip . 

A widely used numerical method for computing the ground state of a single 
component condensate is the imaginary time method followed by an appropriate 
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discretization scheme |15[ [51 13] to evolve the resulted gradient flow equation under 
normalization of the wavefunction, which is mathematically justified by using the 
normalized gradient flow '5', However, it is not obvious that this most popular 
and powerful normalized gradient flow (or imaginary time method) could be directly 
extended to compute the ground state of spin-1 BEC. The reason is that we only have 
two normalization conditions (i.e. the two constraints: conservation of total mass and 
magnetization) which are insufficient to determine the three projection constants for 
the three components of the wavefunction used in the normalization step. In physics 
literatures, the imaginary time method is still applied to compute the ground state of 
spin-1 BEC through the introduction of a random variable to choose the three pro- 
jection parameters in the projection step [501 132] ■ Of course, this is not a determinate 
and efficient way to compute the ground state of spin-1 BEC due to the choice of 
the random variable. Recently, Bao and Wang |5] have proposed a continuous nor- 
malized gradient flow (CNGF) for computing the ground state of spin-1 BEC. The 
CNGF is discretized by Crank-Nicolson finite difference method with a proper and 
very special way to deal with the nonlinear terms and thus the discretization scheme 
can be proved to be mass and magnetization conservative and energy diminishing in 
the discretized level [9]. However, at each time step, a fully nonlinear system must be 
solved which is a little tedious from computational point of view since the CNGF is 
an integral-differential equations (see details in (|A.l[) - (jA.9P ') which involves implicitly 
the Lagrange multipliers in the normalized gradient flow evolution [9] . The aim of this 
paper is to introduce a third normalization condition based on the relation between 
the chemical potentials of spin-1 BEC, in addition to the two existing normalization 
conditions given by the conservation of the total mass and magnetization. Thus we 
can completely determine the three projection constants used in the normalization 
step for the normalized gradient flow. This allows us to develop the most popular and 
powerful normalized gradient flow or imaginary time method to compute the ground 
state of spin-1 BEC. 

The paper is organized as follows. In section [21 the normalized gradient flow 
is constructed by introducing the third projection or normalization condition for 
computing the ground state of spin-1 BEC. In section [31 the backward-forward Eu- 
ler sine-pseudospectral method (BESP) is presented to discretize the normalized 
gradient flow. In section [4l ground states of spin-1 BEC are reported with fer- 
romagnetic/antiferromagnetic interaction and harmonic/optical lattice potential in 
one/three dimensions, respectively. Finally, some conclusions are drawn in section [5l 

2. The normalized gradient flow. In this section, we will construct the nor- 
malized gradient flow for computing the ground state of spin-1 BEC by introducing 
the third normalization condition. 

Various algorithms for computing the minimizer of the nonconvex minimization 
problem (jl.lSp have been studied in literature. For instance, a continuous normal- 
ized gradient flow (CNGF) and its discretization that preserve the total mass and 
magnetization conservation and energy diminishing properties were presented in [Sj. 
Perhaps one of the more popular and efficient techniques for dealing with the nor- 
malization constraints in (|1.16p is through the following construction: choose a time 
step k = At > and denote time steps as tn = n k for n = 0, 1, 2, • • • . To adapt 
an efficient algorithm for the solution of the usual gradient ffow to the minimization 
problem under constraints, it is natural to consider the following splitting (or pro- 
jection) scheme, which was widely used in the physics literature for computing the 
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ground state of BECs: 
'1 



9t(/'i(x,t) 
(2.1) 

9t(/>o(x,t) 

(2.2) 

9t0_i(x,t) 
(2.3) 



y(x) - {(3n +Ps) {\C^1? + I^Z-OP) - {Pn-fisM-, 



-(3s 0-1 00 > 



y(x)-(/3„+/3,)(|0i 



-2/3s 0-1 00 01 



X e 



1. 



l/(x)-(/3„+/3.)(|0-i 



71 > 1, 
- A)|01P 



-/3s 00 01 



X e 



followed by a projection step as 

(2.4) 0i(x,t„) :=0i(x,4) ==c7'i' 0i(x,i-), 

(2.5) 0o(x,i„) := 0o(x,4) = 0o(x,i,7), 

(2.6) 0-i(x,t„) :=0-i(x,4) = a!!i 0-i(x,i-); 

where (j)i{x,t^) = limj_^j± (f)i[:x.,t) {I — —1,0,1) and cr" {I 
constants and they are chosen such that 



-1, 0, 1) are projection 



(2.7) ||<i>(-,t„)|P = \\M;tn)r = 1, Il0l(-,tn)ll' - ll0-l(-,tn)|P = M. 

/=-l 

In fact, the gradient flow (|2.ip - (|2.3p can be viewed as applying the steepest decent 
method to the energy functional E{^) in (|1.14p without constraints, and (|2.4p - (|2.6p 
project the solution back to the unit sphere S in order to satisfy the constraints in 
([TTg)) . In addition, (P?T|) - ((0)l can also be obtained from the coupled GPEs (fO)) - 
(|l.lip by the change of variable t ^ —i t, that is why the algorithm is usually called 
as the imaginary time method in the physics literatures [THl M 13] ■ 
Plugging lOll-dMl) into ^1^, we obtain 



(2.8) 



(2.9) 



i=-i 



Kfll0i(-,<,T)lP-(^-i)'ll0-i(-,^T)lP = M- 



There are three unknowns and only two equations in the above nonlinear system, 
so the solution is undetermined! In order to determine the projection constants a" 
(Z = —1, 0, 1), we need to find an additional equation. Based on the relation between 
the chemical potentials in (|1.23p and the continuous normalized gradient flow proposed 
in [2] for computing the ground state of spin-1 EEC, see details in Appendix A, we 
propose to use the following equation as the third normalization condition 

(2.10) < a-, = {a-f . 

Solving the nonhnear system (|2.8p . (|2.9p and (|2.10p . see details in Appendix B, we 
get explicitly the projection constants as 

(2.11) 

(To - 



l0o(-,tn)|P + ^4(l-M2)||0i(.,t-)||2||0_i(.,i-)||2+M2|l0o(.,i^)r 



1/2 ' 
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(2.12) 



V2Ui{-,tn)\\ ' V2U-l{-,tn)\\ 

From the numerical point of view, the gradient flow (|2.ip - (|2.3p can be solved via 
traditional techniques, and the normalization of the gradient flow is simply achieved 
by a projection at the end of each time step. 

3. Backward-forward Euler sine-pseudospectral method. In this section, 
we will present the backward-forward Euler sine-pseudospectral method (BESP) to 
discretizc the normalized gradient flow ((2TT|l - ((2?3)) . (IT4| - (p?6l) and ((2lT|l - ((27T2)) . 

Due to the trapping potential V^(x) given by (|1.4p . the solution <i>(x, t) decays to 
zero exponentially fast when |x| — > oo. Thus in practical computation, we truncate 
the problem into a bounded computational domain fix (chosen as an interval (a, b) in 
ID, a rectangle (a, 6) x (c, d) in 2D, and a box (a, 6) x (c, d) x (e, /) in 3D, with |a|, 
|c|, |e|, 6, d and / sufficiently large) with homogeneous Dirichlet boundary conditions. 

For simplicity of notation we introduce the method for the case of one spatial di- 
mension {d — 1) defined over the interval (a, h) with homogeneous Dirichlet boundary 
conditions. Generalization to higher dimension are straightforward for tensor product 
grids, and the results remain valid without modifications. For d = 1, we choose the 
spatial mesh size h = Ax > with h — {b — a)/M for M an even positive integer, 
and let the grid points be 

xi := a + j h, j = 0, 1, • • • ,M. 

Let <I>^ = ,,)^ be the approximation of $(a;j,i„) = {(j)i{xj,tn), (l)o{xj,tn), 

0„i (xj, t„))"^ and $" be the solution vector with component In the discretization, 
we use sine-pseudospectral method for spatial derivatives and backward/forward Euler 
scheme for linear/nonlinear terms in time discretization. The gradient flow (|2.ip - (|2.3p 
is discretized, for j = 1, 2, . . . , Af — 1 and n > 1, as 

(3.1) = 2^LriU=., - + G^-\ 

/A* j^n— 1 ^ 

(3.2) ^''^ ^J"'-'- ^ -i?L'/'SU=.. - o^o^o, + Go7,\ , 



where 



G^^ji = [a^ - V{x,) - +/3.) (l-^-^jT + \4,--'\^) - - A)|</>^I,i/] -/.^J^ 
(3.4) -Ps'^-t, {<PZ'f^ 



(3.5) -2/3^ (l>\y 00^^- ^ ^ 



(3.6) (0^ji)' 
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Here, D^^, a pseudospectral differential operator approximation of d^x-, is defined as 

M-l 



X! fJ.li{U)m sm{fim{xj - a)), j = f,2, • • • ,M- f, 



where {U)m (m = 1,2, ■ ■ ■ , M — 1), the sine transform coefficients of the vector U 
{Uq, Ui, ■ ■ ■ , Um)^ satisfying C/q = Um = 0, are defined as 



A/-1 



M 



Uj sm.{^m{xj - a)), m = f,2,- • • ,M- 1; 



and ai (l = — f , 0, f ) are the stabilization parameters which are chosen in the 'optimal' 
form (such that the time step can be chosen as large as possible) as [3] 



(3.7) 
with 



max I Lmin\ 



{bo 



max I Lmin\ 



„ / Lma 



max I Lmm\ 



-1 ; ' 



^max 


— max 

1<J<JW- 




F(/3„+/3,)(|0^ji 




n + (/3„-/?.)|</'!!lir 


^min 


= min 

l<j<M- 




-(/3„ + /3.)(l<,7'l 






7 max 
^0 


= max 

l<j<M- 


_^ [V{x,) ' 


f (/3„+/3,)(|^^-i 






imin 


= min 

l<j<M- 




-(/3„+/3.) (105^-^1 






^max 


= max 

l<j<M- 


^ [V{x,) ' 


f (/3„+/3.)(|0!!7^^,. 


1' + 


|2) + (/3„-A)|^5'jY: 


^min 


= min 

l<j<M- 


JV{x,)^ 


-(/3„+/3.) (10^71 




r) + (/3„-/3,)|^5'-Y] 



The homogeneous Dirichlet boundary conditions are discretized as 

(3.8) 01 4>Im = 00,0 = <^0,M = 4>-UQ = </'-13/ = 0. 

The projection step (|2.4p - (|2.4p is discretized, for < j < M and n> 1, as 
(3.9) 



'1 V^lj: 



where 
(3.10) 



^0 V^Oji 



'-1 



mV + V4(l-M2)||0*p||01^||2 + M2||0*||4 



1/2 ' 



(3.11) 
with 



V2m\\ 



M-l 



Vl-Af-a^||<^Sl| 

V2UU\\ 
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The initial data (|A.10[) is discretized as 



</.?,^- =0i(x,,O), J =0,1,2,... ,M, ; = -1,0,1. 

The linear system p.ip - p.3p can be solved very efficiently by using the fast sine 
transform. In fact, take discrete sine transform at both sides, we get 



(3.12) 
(3.13) 
(3.14) 



1 

At 
1 

At 
1 

At 



{(i>o)r. 



Mm +"1 



'>l)'m + (G" ^)„i, 



i^o)m + (Go ^)r. 



1 < TO < M, 





'1 2 







_l)m + (G™!^), 



Solve the above system in the phase space, we obtain 

1 



(3.15) 
(3.16) 
(3.17) 



1 



At [ai 
1 



M^/2] 



1 + A< [a, + til/2] 



1 + At [a, + fil/2] 



(Gr^)„ 



1 < TO < M, 



Remark 3.1. The gradient flow (K7\)-(K^ can also be discretized by using 
the backward Euler finite difference method proposed in |^ or the backward Euler 
sine-pseudospectral method proposed in ^ for computing the ground state of one- 
component BEC. 

4. Numerical Results. In this section, we ffist show that the ground states 
computed by our new numerical method are independent of the choice of the initial 
data in (jA.lOp and verify numerically the energy diminishing property of the method. 
Finally, we apply the method to compute the ground state of spin-1 BEC with different 
interactions and trapping potentials. In our computations, the ground state is reached 
by using the numerical method ([331) -([SJl), ((3:g|) - ((3lI1) when ||$'^'+^ - < e := 
10~^. In addition, in the ground state of spin-1 BEC, we have M ^ —M (fii ^ 

thus we only present results for < M < 1. 

4.1. Choice of initial data. In our tests, two typical physical experiments are 
considered: 

• Case I. With ferromagnetic interaction, e.g. *^Rb confined in a cigar-shaped 
trapping potential with parameters: to = 1.443 x 10~^^[kg], oq — 5.387[nm], 
a2 = 5.313[nm], cj^ = 27r x 20[Hz], = = 27r x 400[Hz]. This sug- 
gests us to use dimensionless quantities in ()1.9p - ()l.lip for our computa- 
= 1, V{x) = 2:2/2, /3„ w ^"^""at!"''^ = 0.08857V and 

47r(a2-ao)jv ^ -0.000417V with N thc total number of atoms in the 

condensate and the dimensionless length unit Og = ^JhjmLiJx = 2.4116 x 10~^ 
[m] and time unit tg = l/cOx = 0.007958 [s]. 

Case II. With antifcrromagnetic interaction, e.g. ^^Na. confined in a cigar- 
shaped trapping potential with parameters: to = 3.816 x lO^^ejj^g]^ = 
2.646[nm], 02 = 2.911[nm], lu^ ^ 2t: x 20[Hz], = = 27r x 400[Hz]. 
Again, this suggests us to use the following dimensionless quantities in our 
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computations: d = 1, V{x) = l3„ w 0.0241iV and /J^ w 0.000757V 

with the dimensionless length unit = 4.6896 x 10^^ [ni] and time unit 
ts = 0.007958[s]. 

We first test that the converged solution is independent of different choices of 
the initial data in (|A.10[) and energy diminishing property of the normalized gradient 
flow. In order to do so, we choose the initial data in (jA.10[) as 

• Gaussian profiles satisfying the constraints in (|1.16p initially, i.e. 



(4.1) 
(4.2) 
(4.3) 



0i(x,O) 



rl/4 



ri/4' 



-OO < X < oo, 



^-iix,0) 



ri/4 



where k is a constant satisfying < k < 1 — |Af| 
Unnormalized Gaussian profiles, i.e. 



(4.4) 



0i(a;,O) = Mx,0) = <^-i(a;,0) = e-"'/^, 



oo < a; < oo. 



50 100 150 



50 100 
t 



50 100 



Fig. 4.1. Time evolution of Ni = ||(/>i(-,t)p ('left'), No = \\4io{-,t)\\^ ('middle') and N-i = 
||</>_i(-,i)|P ('right') by our method (E4\l- (k^} for ^"^ Rb in Case I with M = 0.5 and N = 10* to 
analyze the convergence of different initial data in [4-4^ (solid line) and | f4.-?p - |?~3| ) voith k = 0.1 
(dotted line), k = 0.2 (dash-dot line) and k = 0.4 (dashed line), respectively. 

We solve the problem (I1.15P by our method on [—16, 16] with time step At = 0.005 
and mesh size h = 1/64 for different values of k in (|4.ip - (|4.3p . Figure [01 plots time 
evolution of Ni{t) ||0i(-,t)||2 [l = 1,0,-1) for ^'^Rb in Case I with M = 0.5 and 
N = 10"* for different choices of the initial data in (|4.4|1 and (|4.1|l - (|4.3p . and Figure 
14.21 shows similar results for ^^Na in Case II. In addition, Figure 14.31 depicts time 
evolution of the energy for the two cases with M — 0.5 and — IQ* for different 
choices of the initial data in (14. 4p . 

From Figs. 14.11 and 14. 2[ we can see that the converged ground states are in- 
dependent of the choice of initial data. In fact, based on our extensive numerical 
experiments on other types of initial data (not shown here for brevity) , our numerical 
method always gives the ground state if all the three components in the initial data 
are chosen as nonnegative functions. In addition. Fig. 14.31 demonstrates the energy 
diminishing property of the normalized gradient flow and its full discretization when 



Computing ground states of spin-1 BEC by normalized gradient flow 



11 




02468 10 02468 10 '02468 10 



Fig. 4.2. Time evolution of Ni = ||<^>i (■, t) |p ('left'), Nq = ||0o(-,t)lP ('middle') and Af_i = 
||</>_l(-,t)|p ('right') by our method (E^p-El) for ^^Na in Case II with M = 0.5, and N = 10'' to 
analyze the convergence of different initial data in \4-4^ (solid line) and with k = 0.1 

(dotted line), k = 0.2 (dash-dot line) and k = 0.4 (dashed line), respectively. 




Fig. 4.3. Time evolution of the energy by our method \2.4}^ - li2.6\) with M = 0.5, and N = 10* 
for a) ^"^ Rb in case I; and b) Na in case II with different initial data in (solid line) and 

^T7p-{73^ with K = 0.1 (dotted line), k = 0.2 (dash-dot line) and k = 0.4 (dashed line), respectively. 

time step Ai is small. Based on our numerical experiments, for < M < 1, we 
suggest the initial data in (jA.lOp be chosen as: i) with ferromagnetic interaction, i.e. 

/?. < 



<^i(x) = -VTT3M07(x), 0o(x) = ^^— 07(x), 0i(x) = -VT^07(x); 
and ii) with antiferromagnetic interaction, i.e. /3s > 

</.i(x) = yi±^07(x), </.o(x)=0, 0i(x) = yi^0^P(x); 

where </'gP(x) can be chosen as the approximate ground state solution of single com- 
ponent BEC, e.g. the harmonic oscillator approximation when /3„ small and the 
Thomas- Fermi approximation when /3„ ^ 1 [51 [3 [7] . Based on these choices of initial 
data, we report the ground states computed by our numerical method. 

Figure 14.41 shows the ground state solutions of ®^Rb in Case I with N — 10^ 
for different magnetization M and Table 14.11 lists the corresponding ground state 
energies and their Lagrange multipliers (see their detailed formulation in Appendix 
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Fig. 4.4. Wave functions of the ground state, i.e. <7!>i(a;) (dashed line), (poi^) (solid line) 
and tp-ii^) (dotted line), of ^"^ Rb in Case I with fixed number of particles N = 10^ for different 
magnetization M = 0, 0.2, 0.5, 0.9. 

Similarly, Figure 14.61 shows the ground state solutions of ^■^Na in Case II with 
N = 10* for different magnetization M and Table l472l lists the corresponding ground 
state energies and their Lagrange multipliers. In addition. Figure [577] shows similar 
ground state solutions with M — 0.5 for different particle number N . 

Figure 14.81 plots the mass of the three components in the spin-1 BEC ground 
states with N = 10* for different magnetization M, and Figure depicts the energy 
and chemical potentials with M = 0.5 for different particle number N. 

From Figs. I4.4ll4.6l as well as Tabs. 14.1114.21 we can draw the following conclusions: 
(i) For ferromagnetic interaction in the spin-1 BEC, i.e. f3s < 0, the three components 
in the ground state solutions are all positive functions (c.f. Figs. 14.41 and l4.5p ; while 
for antiferromagnetic interaction, i.e. /3s > 0, (f>i and cjj-i are positive functions and 
00 = (c.f. Figs. 14.61 and l4.7p . (ii) For ferromagnetic interaction in the spin-1 BEC, 
i.e. Ps < 0, for fixed number of particles N in the condensate, when the magnetization 
M increases from to 1, the mass A'^i increases from 0.25 to 1, the mass A^_i decreases 
from 0.25 to and the mass Nq decreases from 0.5 to (c.f. Fig. 14. 9h ): while for 
antiferromagnetic interaction, i.e. f3s > 0, Ni increases from 0.5 to 1, A^_i decreases 
from 0.5 to and A^o = (c.f. Fig. l4.9b V (iii) For ferromagnetic interaction in 
the spin-1 BEC, i.e. /3s < 0, for fixed number of particles N in the condensate. 
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M 




E 






A(xlO- 





36 


1365 


60 


2139 





0.1 


36 


1365 


60 


2139 


1.574 


0.2 


36 


1365 


60 


2139 


1.621 


0.3 


36 


1365 


60 


2139 


1.702 


0.4 


36 


1365 


60 


2139 


1.827 


0.5 


36 


1365 


60 


2139 


2.014 


0.6 


36 


1365 


60 


2139 


2.218 


0.7 


36 


1365 


60 


2139 


2.062 


0.8 


36 


1365 


60 


2139 


2.081 


0.9 


36 


1365 


60 


2139 


2.521 



Table 4.1 

Ground state energy E and their chemical potentials fi and X for ^''R6 in Case I with N = 10* 
for different magnetization M. 



N=10' 



N=10'' 
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12 -12 
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Fig. 4.5. Wave functions of the ground state, i.e. </)i(x) (dashed line), 4>o(x) (solid line) and 
if>—-i{x) (dotted line), of^'^Rb in Case I with magnetization M = 0.5 for different number of particles 
N. 



the energy and chemical potentials are almost independent of the magnetization (c.f. 
Tab. 14. 1|) : while for antifcrromagnetic interaction, i.e. f3s > 0, when the magnetization 
M increases from to 1, the energy E increases and the main chemical potential fj, 
decreases and the second chemical potential A increases (c.f. Tab. I4.2p . In both cases, 
for fixed magnetization M , when the number of particles N increases, the energy and 
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M=0 



M=0.2 




Fig. 4.6. Wave functions of the ground state, i.e. 4>i{x) (dashed line), (poi^) (solid line) 
and (j)—i{x) (dotted line), of Na in Case II with fixed number of particles N = lO** for different 
magnetization M = 0,0.2,0.5,0.9. 



chemical potentials increase (c.f. Fig. 14. 8p . These observations agree with those 
obtained in ^ and f3D^ by different numerical methods. 

4.2. Application in ID with optical lattice potential. In this subsection, 
our method is applied to compute the ground state of spin-1 BEC in one dimension 
(ID) with an optical lattice potential. Again, two different interaction are considered: 

• Case I. For *^Rb with dimensionless quantities in (|1.9p - (|l.lip used as: d = 1, 
V{x) = + 25 sin^ (^) , ^„ = 0.0885A and = -0.00041 A, with A the 
total number of atoms in the condensate and the dimensionless length unit 
as = 2.4116 X 10-^ [m] and time unit ts = 0.007958 [s]. 

• Case II. For ^'^Na with dimensionless quantities in ()1.9p - ()l.lip used as: d = 1, 
V{x) = + 25sin2 {^), /3„ = 0.0241 A and = 0.00075A, with A the 
total number of atoms in the condensate and the dimensionless length unit 
Us = 4.6896 X 10-*' [m] and time unit = 0.007958[s]. 

Figure l4. 101 shows the ground state solutions of ^^Rb in Case I with A — lO"' for 
different magnetization M and Table lists the corresponding ground state energies 
and their Lagrange multipliers. Figure l4TTT] and Table [44] show similar results for ^'^Na 
in Case II. 

From Figs. l4.10Ej^?TT] and Tabs. l4.3Bj4.41 it can be seen that our method can 
be used in computing ground state of spin-1 BEC with general potential. In addition 
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M 


E 




A* 


A 





15.2485 


25, 


,3857 





0.1 


15.2514 


25, 


,3847 


0.0569 


0.2 


15.2599 


25, 


,3815 


0.1142 


0.3 


15.2743 


25, 


,3762 


0.1725 


0.4 


15.2945 


25, 


,3682 


0.2325 


0.5 


15.3209 


25, 


,3572 


0.2950 


0.6 


15.3537 


25, 


,3423 


0.3611 


0.7 


15.3933 


25, 


,3220 


0.4326 


0.8 


15.4405 


25, 


,2939 


0.5121 


0.9 


15.4962 


25, 


,2527 


0.6049 



Table 4.2 

Ground state energy E and their chemical potentials fi and X for Na in Case II with N = 10* 
for different magnetization M . 



ME n A(xlO~^) 






47.6944 


73.0199 





0.1 


47.6944 


73.0199 


0.711 


0.2 


47.6944 


73.0199 


0.788 


0.3 


47.6944 


73.0199 


0.859 


0.4 


47.6944 


73.0199 


0.948 


0.5 


47.6944 


73.0199 


1.072 


0.6 


47.6944 


73.0199 


1.178 


0.7 


47.6944 


73.0199 


1.164 


0.8 


47.6944 


73.0199 


1.200 


0.9 


47.6944 


73.0199 


1.477 



Table 4.3 

Ground state energy E and their chemical potentials and A for Rb in Case I with N = 10* 
for different magnetization M in an optical lattice potential. 



to that, similar conclusions as those in the end of previous subsection can also be 
observed in this case. 

4.3. Applications in 3D with optical lattice potential. In this subsection, 
our method is applied to compute the ground state of spin-1 BEC in three dimensions 
(3D) with an optical lattice potential. Again, two different interaction are considered: 

• Case I. For ^^Rb with dimensionless quantities in (|1.9p - (|l.lip used as: d = 
3, V{x) = l[x^+y^ + z^) + 100 [sin^ {^) + sin^ (^) + sin^ {^)] , /?„ = 
0.0880iV and = -0.00041 A, with N the total number of atoms in the 
condensate and the dimensionless length unit Os = ^Jh/mojx = 7.6262 x 10~^ 
[m] and time imit tg = 1/ujx = 7.9577 x 10^"* [s] (corresponding to physical 
trapping frequencies uj^ = ujy = ujz = 'iir x 200 [Hz]). 

• Case II. For ^^Na with dimensionless quantities in (|1.9p - (|l.lip used as: d = 
3, V{x) = i (a;2 + y2 + + iQO [sin^ (if) + sin^ (^) + sin^ {^)] , /?„ = 
0.0239A and Ps = 0.00075A with A the total number of atoms in the con- 
densate and the dimensionless length unit Ug = 1.4830 x 10^^ [m] and time 
unit tg — 7.9577 x 10^^ [s] (corresponding to physical trapping frequencies 
LUx — LUy ~ LUz ~ 2Tr X 200 [Hz]). 
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Fig. 4.7. Wave functions of the ground state, i.e. 4>i{x) (dashed line), (poi^) (solid line) 
and (j)—i{x) (dotted line), of ^"^ Na in Case II with magnetization M = 0.5 for different number of 
particles N . 



Figure [4. 121 shows the ground state solutions with N = 10^ and M = 0.5 for the 
two cases. 

From Fig. 14.121 we can see that our method can be used to compute the ground 
state of spin-1 BEC in 3D with general trapping potential. 

5. Conclusions. We have proposed an efficient and accurate normalized gra- 
dient flow or imaginary time method to compute the ground state of spin-1 Bose- 
Einstein condensates by introducing a third normalization condition, in addition to 
the conservation of total particle number and the conservation of total magnetization. 
The condition is derived from the relationships between the chemical potentials of the 
three spinor components together with a splitting scheme applied to the continuous 
normalized gradient flows proposed to compute the ground state of spin-1 BEC. The 
backward-forward sine-pseudospectral method is applied to discretize the normalized 
gradient flow for practical computation. The ground state solutions and fraction of 
each component are reported for both ferromagnetic and antiferromagnetic interac- 
tion cases. The energy and chemical potentials of the condensate are also reported. In 
addition, the method may be further extended to other spinor condensate with higher 
degree of freedom as well as spinor condensate in the presence of external magnetic 
field, which will be our future study. 

Finally, based on our extensive numerical experiments and results, we conjecture 
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Fig. 4.9. Energy E and chemical potentials fi and X of spin-1 BEC with fixed magnetization 
M = 0.5 for different number of particles N. a) for in case I; and b) for '^'^ Na in case II. 



that when /3„ > 0, /3„ > \(3s\ and V{x) > satisfying lim^^^^^ V {x) oo, there 
exists minimizer of the nonconvex minimization problem (|1.15|) . In addition, when 
(3s < 0, positive minimizer (the three components are positive function) is unique; 
when (3s > 0, nonnegative minimizer {(f>i and are positive and 0o = 0) is unique. 
Rigorous mathematical justification are on-going. 

Acknowledgment. We acknowledge support from the Ministry of Education of 
Singapore grant No. R-146-000-083-112. 

Appendix A Derivation of the third projection equation (|2.10p 

In order to find the third projection or normalization equation used in the projec- 
tion step of the normalized gradient flow, we first review the continuous normalized 
gradient flow (CNGF) constructed in [3] for computing the ground state of spin-1 
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Fig. 4.10. Wave functions of the ground state, i.e. </)i(x) (dashed line), (poi^) (solid line) 
and tp—ii^) (dotted line), of Rb in Case I with fixed number of particles N = 10^ for different 
magnetization M = 0,0.2,0.5,0.9 in an optical lattice potential. 



BEC in (ITTa : 

9t0i(x,t) = 
(A.l) 



-l3s4>~i(l>l + Mt) + \.i{t)] (1)1, 



dt(f)()i:x-,t) 



(A.2) 



at(/)_i(x,t) 
(A.3) 



- y(x) - (/?„ + + |(/.oP) - (/3n - /3.)|0i|' 

01 -I- [^<i,(t) - A<i.(t)]0_i. 



fi,s>{t) and A$(i) are chosen such that the above CNGF is mass (or normaUzation) and 
magnetization conservative and they are given as [S] 



7V$(i)i?*(t)-Af|(t) 



7V$(t)i?$(i)-M|(i) ' 
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Fig. 4.11. Wave functions of the ground state, i.e. (pi(x) (dashed line), (poi^) (solid line) and 
-i(x) (dotted line), of^^Na in Case 11 with N = 10* for different magnetization M = 0, 0.2, 0.5, 0.9 
an optical lattice potential. 



Fig. 1.12. Contour plots for the wave functions of the ground state, i.e. (pi{x, y, 0) (top row), 
(j)o{x,y,0) (middle row) and ij)-i{x,y,0) (bottom row) with N = 10** and M = 0.5 in an optical 
lattice potential. Left column: for Rb in Case I; and right column: for "^^Na in Case II. 



with 
(A.5) 

(A.6) 

(A.7) 



(A.8) 



M^{t)= I [|</.i(x,t)|2-|0_i(x,i)|2] dx, 

R^{t)= [ [|</)i(x,t)p + |0_i(x,t)|2] rfx, 
Jr'' 

+213, {^-i(t>l^i + <t^-i4>l(t>i) \ dx, 
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37.3775 


0.3458 
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37.3492 


0.4252 
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37.3167 


0.5079 
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25.8696 


37.2305 


0.6920 
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25.9340 


37.2305 


0.6920 



Table 4.4 

Ground state energy E and their chemical potentials fi and X for Na in Case II with N = 10* 
for different magnetization M in an optical lattice potential. 



(A.9) 



+ (/3„ + Ps) |0i|' - + |0oP - 



dx. 



For the above CNGF, for any given initial data 

(A.IO) $(x,0) = (0i(x,O),.^o(x,O),0_i(x,O))^ :=$(°)(x), x e 

satisiying 

(A.ll) N^{t = 0) := N^m ^ 1, AU{t = 0) := 7\f<j,(o) = M, 

it was proven that the total mass and magnetization are conservative and the energy 
is diminishing [9], i.e. 

iV<t,(i) = l, M^{t) = M, £;($(-,t) s)), foranyt>s>0. 

The normalized gradient flow (|2.ip - (|2.6p can be viewed as applying a time-splitting 
scheme to the CNGF (|A.l|) - (|A.3p and the projection step ((T4)l - ((2?6)) is equivalent to 
solving the following nonlinear ordinary differential equations (ODEs): 



(A.12) 
(A.13) 
(A.14) 



9t0i(x,<) = [/i$(t) + A$(t)] 
dt(l)-i{x,t) = [/i$(t) - A$(t)] 



n> 1. 



The solution of the above ODEs can be expressed as 



(A.15) 
(A.16) 
(A.17) 



i'i(x,t„) = exp 



[^$(r) + A$(r)] dr] 0i(x, i„„i), 



6o(x,t„) = exp 



/i$(t) dr 0o(x, 



tr^-l 



-i(x,t„) = exp 



[^$(r) - A$(t)] dr 0_i(x, i„_i). 
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This solution suggests the following relation between the coefficients 



exp ^ J [^i^{T) + A$(t)] drj exp ^ J [n<i,{T) - A$(t)] dr 

(A. 18) = exp 2^^{t) dr^ — exp ^J■■^>iT) dr 

This immediately suggests us to propose the third normalization equation (j2.10p to 
determine the projection parameters. In fact, equation (|2.10p can be also obtained 
from the relation between the chemical potentials in (jl.23p by physical intuitions. 

Appendix B Derivation of the projection parameters in (|2.11|) - (|2.12|) 

Summing (|2TT|) and ([2?T2l) . we get 

(B.l) 2(arf||0i(,t-)|p.= l + M-K)2||^o(-,^T)|p. 

This immediately implies 



Jl + M-iamoi;tn: 

(B.2) < = 



V2\\M;tn)\\ 
Subtracting ([^?T^ from ^J^, we obtain 

(B.3) 2(a«i)2||0_i(.,t-)||2 = i-Af-Kf|l</.o(-,^T)f. 

Again, this immediately implies 



(B.4) cr"i = — 



V2U-l{;tn)\\ 

Multiplying (|B.2|) and (|B.4p and noticing (PTTUl) . we get 

(B.5) =.4||0_i(,t-)f ||0i(,t-)|p(ao")^ 

Simplifying the above equation, we obtain 

[ii^o(-,t;^)r-4ii0„i(.,tjini0i(-,i-)ir] {a^r-2\\M;t-)\\'ia^^r 

(B.6) +(1 - M^) = 0. 

Solving the above equation and noticing (ctq)^ II</'o('5 )|P < (1 ^ ^l'^), we get 

\\M;t-W - V'4(l-^^')ll<^l(-'^«)IPII'/'-l(-'*")P+^^'ll<^0(-,tn)|| 



(B.7) = 



||0o(-,i«)r-4||0-i(-,i^)P ||0i(-,t;:)P 

1 - M2 



ll0o(-,irT)||2 + Y/4(l-M2)||0i(.,i-)||2||0_i(.,t-)||2 + M2||0o(-,t^; 
Thus immediately implies the solution in (j2.1ip . 
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Appendix C Computing the chemical potentials /i and A 

After we get the ground state $ numerically, the energy of the ground state 
can be computed from the discretization of (|1.14p immediately. In order to compute 
the chemical potentials numerically, different formulations can be applied. Here we 
propose one of the most reliable way to compute them. Multiplying both sides of 
(jl.lSp by 01 and integrate over R'^, we get 



(C.l) 



(m + A)||0i 



61 Hi(l)i dx 



L,Hl0i). 



Similarly, take the same procedure to (|1.19p and (|1.20|) by multiplying (f)Q and 
respectively, we obtain 



4>o Hq4>o dx := (00, -ffo^o), 



(C.2) 

(C.3) (a*- A)||0_i||2 = / 0_ii/_i0_idx:=(0_i,i7_i0_i). 



Summing (jC.ip . (|C.2p and (jCSp . noticing that the ground state $ satisfying the 
constraints (|1.16p . we get 

(C.4) /i + Af A= (0l,7^l0l) + (0o,i^o0o) + (0-l,i^-l0-l)■ 

Subtracting (|C.3P from (jCip . we get 

(C.5) Af (Il0iir + II0-1II') A = (0i,i?i0i) - (0-i,i?-i0-i). 



Solving the linear system (IC.4P and (IC.5P for the chemical potentials /i and A as 
unknowns and integrating by parts to the right hand sides, we have 



(C.6) M 
where 



{\\M' + U^ir)D{^)-M F{^) 



Af2 



A = 



F($) - M £)($) 

\\M' + U-iP-M^' 



Pn\<l>0\ 



+ {(3n+f3s) [|0i|" + + 2|0oP (|0l|' + 10- 



(C.7) 



+2A 



i0o'/'i + ^-i0o'/'i) f dx. 



jRd I ^ 



p-|V0_ip)+nx)(|0ip-|0-iP) 



(C.8) 



(/3„ + /?,) |0i|4 - l^^il'i + |0oP (|0iP - 



dx. 



Thus the chemical potentials /i and A can be computed numerically from the dis- 
cretization of dCJl, (|CJl) and JCJl. 
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